StatisticalCode_ModellingMasting

Libraries used

library(tidyverse)
library(hablar)
library(patchwork)
library(plotly)
library(ggparallel)
library(ggforce)
library(lubridate)
library(fmsb)
library(rstatix)

#functions
sd.prop <- function(x, n){
  sqrt((x*(1-x))/n)
}

se <- function(x){
  sd(x)/sqrt(14*729) #n=No of years*number of parameter combinations for 2004-2018
}

shap_test <- function(df){
  shapiro.test(df$med_DBH)
}

shap_test2 <- function(df){
  shapiro.test(df$aver)
}

wil_test <- function(df){
  wilcox.test(med_DBH~PFT, df, alternative='less')
}

wil_testG <- function(df){
  wilcox.test(aver~PFT, df, alternative='greater')
}

#Not significantly different
wil_testNS <- function(df){
  wilcox.test(aver~PFT, df)
} 

#dropping year function
drop_y <- function(df){
  select(df,-Year)
}
#Kappa function
fun_Kap2 <- function(x){
  Kappa.test(x)
}
#getting Kappa values
est_Kap <- function(x){
  x$Result$estimate
}
#getting p values
p_Kap <- function(x){
  x$Result$p.value
}
#effect size
wil_test_eff <- function(df){
  df %>% pivot_longer(cols=c(freq,freqW), names_to='type', values_to='No') %>% 
    wilcox_effsize (No~type)
} 

Analysis and statistical code is shown for RB+Pollen (named as Polim_)

#Field data/Obs Monitoring1000 (70m*150m)
uryu<- read.csv('uryu.csv',row.names = 1,
                colClasses = c(Family='factor', 
                               Latin_name='factor',
                               PFT='factor'))
#Monitoring1000 for 1 hectare, biomass in kg
uryuYear <- read.csv('uryuYear.csv',row.names = 1)

#Seed
#1983-
Yacorns <- read.csv('Yacorns.csv', row.names = 1)
#2004
Yseed <- read.csv('Yseed.csv', row.names = 1,
                  colClasses = c(Latin_name='factor',
                                 mastingW='factor',
                                 mastingN='factor'))

#Simulation data
#QuercusNo and FloweringNo at the time of flowering, Masting and Masting_n at end of year
#MassRatio and FloweringRatio at time of flowering
Polim <- read.csv('Polim.csv', row.names = 1)

Table 1 statistics

#Field data stats
uryuYear %>% 
  summarise(meanNo=mean(TreeNo), sdNo=sd(TreeNo),
            meanNoProp=mean(TreeProp), seNoProp=mean(sdNoProp),
            meanBiom=mean(Biomass), sdBiom=sd(Biomass),
            meanBiomProp=mean(BiomassProp), seBiomProp=mean(sdBiomProp))

#Simulation between 2005-2018
#RB+Polim
#tree number
Polim %>% filter (Year>2004) %>% 
  group_by(Rc, Li, Bu, Year) %>% #will calculate mean from 5 runs for every year
  summarise(Tree=mean(TreeNo), Percent=mean(NoRatio)) %>% ungroup() %>% 
  mutate(sePerc=sd.prop(x=Percent, n=Tree)) %>% 
  summarise(mean=mean(Tree), se=se(Tree), meanPerc=mean(Percent), sePercM=mean(sePerc))

#Biomass
Polim %>% filter (Year>2004) %>% group_by(Rc, Li, Bu, Year) %>% 
  summarise(biomass=mean(WoodyPftMass), Quercus=mean(Masting), Percent=Quercus/biomass,
            Tree=mean(TreeNo), sePerc=sd.prop(x=Percent, n=Tree)) %>%
  ungroup() %>% 
  summarise(mean=2*mean(biomass), se=2*se(biomass), meanQ=2*mean(Quercus), seQ=2*se(Quercus),
            PercentM=mean(Percent), PercentSE=mean(sePerc))

Observation-Simulation comparison

#Obs
uryu %>% group_by(PFT) %>%
  summarise(biomassPft=mean(biomass),
            DBHmin=range(dbh)[1],
            DBHmax=range(dbh)[2],
            DBHmed=median(dbh),
            DBHav=mean(dbh)) #Masting is 642 kg

#normality of biomass and dbh
uryu2 <- uryu %>%
  filter(PFT %in% c('Masting', 'TeBS')) %>% 
  group_by(Year,PFT) %>% summarise(MeanB=mean(biomass), MedBH=median(dbh)) 
tapply(uryu2$MeanB,uryu2$PFT, shapiro.test)#both normal

## $BoBS
## NULL
## 
## $BoNE
## NULL
## 
## $Masting
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.92512, p-value = 0.2604
## 
## 
## $TeBS
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.98663, p-value = 0.997

tapply(uryu2$MedBH,uryu2$PFT, shapiro.test)#both normal

## $BoBS
## NULL
## 
## $BoNE
## NULL
## 
## $Masting
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.93776, p-value = 0.3904
## 
## 
## $TeBS
## 
##  Shapiro-Wilk normality test
## 
## data:  X[[i]]
## W = 0.96244, p-value = 0.763

#t-test
t.test(MeanB~PFT, data=uryu2) #mean biomass

## 
##  Welch Two Sample t-test
## 
## data:  MeanB by PFT
## t = 40.96, df = 13.091, p-value = 3.297e-15
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  508.4844 565.0672
## sample estimates:
## mean in group Masting    mean in group TeBS 
##              644.1987              107.4229

t.test(MedBH~PFT, data=uryu2) #median DBH

## 
##  Welch Two Sample t-test
## 
## data:  MedBH by PFT
## t = 5.7635, df = 21.297, p-value = 9.6e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.674635 1.435305
## sample estimates:
## mean in group Masting    mean in group TeBS 
##              11.11015              10.05518

#Sim
#RB+Pollen-median DBH
trial2 <- Polim %>% select(Rc:run, TeNE_medDBH:BoBS_medDBH) %>% 
  pivot_longer(cols=TeNE_medDBH:BoBS_medDBH, names_to='PFT', values_to='median') %>% 
  group_by(Year, PFT, Rc, Li, Bu) %>% 
  summarise(med_DBH=mean(median)) #eliminating runs

#normality
#nesting
trial2  %>%  mutate(med_DBH=med_DBH+rnorm(n(),mean = 0, sd=0.0001))%>% #correcting for all equal problem
  group_by(PFT, Rc, Li, Bu) %>% 
  nest() %>% 
  mutate(shapTest=map(data, shap_test), glance=map(shapTest, broom::glance)) %>% 
  unnest(glance) %>% filter(p.value<0.05) %>%  
  group_by(PFT) %>% count() #Mostly non-normal

#Wilcox test
trial_Wil <- Polim %>% select(Rc:run, TeNE_medDBH:BoBS_medDBH) %>% 
  pivot_longer(cols=TeNE_medDBH:BoBS_medDBH, names_to='PFT', values_to='median') %>%
  filter(PFT %in% c('TeBS_medDBH','Masting_medDBH')) %>% 
  mutate(PFT=fct_relevel(PFT, 'TeBS_medDBH')) %>% 
  group_by(Year, PFT, Rc, Li, Bu) %>% summarise(med_DBH=mean(median))
trial_Wil %>% group_by(Rc, Li, Bu) %>% 
  nest() %>% 
  mutate(WilTest=map(data, wil_test), glance=map(WilTest, broom::glance)) %>% 
  unnest(glance) %>% filter(p.value < 0.05) #Masting PFT never higher

Mean biomass and selection of viable parameter combinations

#RB+Pollen Mean Biomass
Polim_stat <- Polim %>% mutate(TeNE=TeNE/TeNE_n,
                               TeBS=TeBS/TeBS_n,
                               BoNE=BoNE/BoNE_n,
                               BoNS=BoNS/BoNS_n,
                               BoBS=BoBS/BoBS_n,
                               Masting=Masting/Masting_n) %>% 
  select(Rc:run, TeNE:Masting) %>% 
  pivot_longer(cols=TeNE:Masting, names_to='PFT', values_to='aver') %>% 
  filter(PFT %in% c('TeBS','Masting'))
#nesting
Polim_stat  %>% group_by(PFT, Rc, Li, Bu) %>%  #taking together 5 runs Masting or TeBS
  nest() %>% 
  mutate(shapTest=map(data, shap_test2), glance=map(shapTest, broom::glance)) %>% 
  unnest(glance) %>% filter(p.value<0.05) %>%  
  group_by(PFT) %>% count() #Mostly non-normal

#Getting parameter combinations, Masting S >
Polim_par_g <- Polim_stat %>% group_by(Rc, Li, Bu, Year, PFT) %>% summarise(aver=mean(aver)) %>% 
  group_by(Rc, Li, Bu) %>% 
  nest() %>% 
  mutate(WilTest=map(data, wil_testG), glance=map(WilTest, broom::glance)) %>% 
  unnest(glance) %>% filter(p.value<0.05)  %>% 
  unite('parameter', c('Rc', 'Li','Bu'), sep='_', remove = F) %>% .$parameter
#NS
Polim_par_e <- Polim_stat %>% group_by(Rc, Li, Bu, Year, PFT) %>% summarise(aver=mean(aver)) %>% 
  group_by(Rc, Li, Bu) %>% 
  nest() %>% 
  mutate(WilTest=map(data, wil_testNS), glance=map(WilTest, broom::glance)) %>% 
  unnest(glance) %>% filter(p.value>=0.05) %>% 
  unite('parameter', c('Rc', 'Li','Bu'), sep='_', remove = F) %>% .$parameter

#parameters unite
Polim_par <- union(Polim_par_e, Polim_par_g)

#filtered dataset
Polim_f <- Polim %>% 
  unite('parameter', c('Rc', 'Li','Bu'), sep='_', remove = F) %>% 
  filter(parameter %in% Polim_par) 

Flowering CV

Polim_f %>% filter(Year>2004) %>% 
  group_by(Rc, Li, Bu, run) %>% summarise(FlowerCV=sd(FloweringRatio)/mean(FloweringRatio)) %>% 
  group_by(Rc, Li, Bu) %>% summarise(meanFlowerCV=mean(FlowerCV), seFlowerCV=sd(FlowerCV)/sqrt(5)) %>%
  arrange(desc(meanFlowerCV)) 

Figure 3

#Uryu
e <- uryu %>% group_by(Year, PFT) %>% 
  summarise(biomass=sum(biomass), count=n()) %>% #for every Year, every family average
  mutate(aver=biomass/count) %>% 
  ggplot(aes(x=PFT, y=aver))+
  geom_boxplot()+
  labs(title='(e)', y='Mean biomass/tree (kg dry mass)')+
  theme_bw()

#RB+Pollen
c <- Polim %>% mutate(TeNE=TeNE/TeNE_n,
                         TeBS=TeBS/TeBS_n,
                         BoNE=BoNE/BoNE_n,
                         BoNS=BoNS/BoNS_n,
                         BoBS=BoBS/BoBS_n,
                         Masting=Masting/Masting_n) %>% 
  select(Rc:run, TeNE:Masting) %>% 
  pivot_longer(cols=TeNE:Masting, names_to='PFT', values_to='aver') %>%
  filter(Year>2004) %>% group_by(Year, PFT, Rc, Li, Bu) %>% 
  summarise(aver=mean(aver)) %>% 
  ggplot(aes(x=PFT, y=aver*1000*2))+
  geom_boxplot()+
  labs(title='(c)', y='Mean biomass/tree (kg dry mass)')+
  theme_bw()
c+e

Figure 4

#RB+Pollen
Polim_mast_meanBiom <- Polim %>% mutate(TeNE=TeNE/TeNE_n,
                                        TeBS=TeBS/TeBS_n,
                                        BoNE=BoNE/BoNE_n,
                                        BoNS=BoNS/BoNS_n,
                                        BoBS=BoBS/BoBS_n,
                                        Masting=Masting/Masting_n) %>% 
  select(Rc:run, TeNE:Masting) %>% 
  pivot_longer(cols=TeNE:Masting, names_to='PFT', values_to='aver') %>%
  filter(Year>2004) %>% group_by(Year, PFT, Rc, Li, Bu) %>% 
  summarise(aver=mean(aver)) %>% filter(PFT=='Masting')
trial <- Polim_mast_meanBiom %>% filter(Rc==5, Year==2010) %>% ungroup() %>% 
  select(Bu, Li, aver) %>% mutate(aver=aver*1000*2) %>% 
  convert(num(Bu, Li)) %>% mutate(Bu=Bu*100, Li=Li*100) 
trial2 <- xtabs(aver~Li+Bu, data=trial) #9*9
trial3 <- as.numeric(trial2)
trial4 <- matrix(trial3, nrow = 9, ncol=9)  
plot_ly(x=as.numeric(colnames(trial2)),
             y=as.numeric(rownames(trial2)),
             z=trial4, type = 'surface', colorscale = list(c(0, 1), c('white', "blue")),
             cmin=0, cmax=220,
             contours=list(
               z=list(show = TRUE, start = 0.0, end = 220, size = 6, color = 'black')
             )) %>% 
  layout(title=list(text=paste0('(c)'), y=0.97, x=0.2),font=list(size=14),
         scene=list(aspectratio=list(x=1,y=1,z=1),
                    xaxis=list(title='<i>B<sub>u</sub></i> (%)'),
                    yaxis=list(title='<i>L<sub>T</sub></i> (%)'),
                    zaxis=list(title='Mean biomass (kg)',range=c(0,220)),
                    camera = list(eye = list(x = 1.3, y = 1.3, z = 1.3),
                                  center = list(x=0, y=0, z=-0.2)))) %>% hide_colorbar()

Figure 5

#other filtered simulations
Basic_f <- read.csv('Basic_f.csv',row.names = 1)
BasicW_f <- read.csv('BasicW_f.csv',row.names = 1)
PolimW_f <- read.csv('PolimW_f.csv',row.names = 1)

#parameters filtered for each masting function
Basic_f <- Basic_f %>%  mutate(Type='RB')
BasicW_f <- BasicW_f %>%  mutate(Type='RB-W')
Polim_f <- Polim_f %>%  mutate(Type='RBP')
PolimW_f <- PolimW_f %>%  mutate(Type='RBP-W')

#Combining
combined <- bind_rows(Basic_f, BasicW_f, Polim_f, PolimW_f)
combined <- combined %>% select(Rc:Bu, Type) %>% convert(fct(Type)) %>% 
  mutate(Type=fct_relevel(Type, 'RB', 'RBP', 'RB-W'))
combined2 <- combined %>% convert(num(Rc:Bu)) %>% 
  mutate(Li=Li*100, Bu=Bu*100) %>% 
  group_by(Rc, Li, Bu, Type) %>% count() 

combined3 <- gather_set_data(combined2,1:3)

combined3 %>% mutate(Type=fct_recode(Type, '(a)'='RB',
                                     '(b)'='RB-W',
                                     '(c)'='RBP',
                                     '(d)'='RBP-W')) %>%
  ggplot(aes(x, id = id, split = y, value = n)) +
  geom_parallel_sets( alpha = 0.3, axis.width = 0.1) +
  geom_parallel_sets_axes(axis.width = 0.3, fill='white', color='grey') +
  geom_parallel_sets_labels(colour = 'black' , angle = 0)+
  facet_wrap(~Type)+
  theme_minimal()+
  theme(axis.text.y=element_blank(),
        axis.ticks = element_blank(),
        axis.text.x = element_text(size=13, face='bold'),
        panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        strip.text.x = element_text(size=12, hjust = 0.1))+
  scale_x_discrete (name='Parameters',
                    labels = c(expression(paste(italic('B'[u]),' (%)')),
                               expression(paste(italic('L'[T]),' (%)')),
                               expression(paste(italic('R'[c])))))

Figure 6

# RB+Pollen
trialPolBu3 <- Polim_f %>% filter(Rc==3, Bu==0.03, Year>1983) %>% 
  group_by(Year, Li) %>% summarise(FloweringRatio=mean(FloweringRatio)) %>% 
  ungroup() %>% 
  select(Year, Li, FloweringRatio) %>% 
  convert(num(Li)) %>% mutate(Li=Li*100) 
trialPolBu4 <- Polim_f %>% filter(Rc==3, Bu==0.04, Year>1983) %>% 
  group_by(Year, Li) %>% summarise(FloweringRatio=mean(FloweringRatio)) %>% 
  ungroup() %>% 
  select(Year, Li, FloweringRatio) %>% 
  convert(num(Li)) %>% mutate(Li=Li*100) 

plot_ly(data=trialPolBu3[which(trialPolBu3$Year>2004),], x = ~Li, y = ~Year, z = ~FloweringRatio, split = ~Li, type = 'scatter3d', 
        mode = 'lines',
        line = list(color='black',width = 2)) %>%  #works really well, different color for lines
  add_trace(data = trialPolBu4[which(trialPolBu4$Year>2004),],
            x = ~Li, y = ~Year, z = ~FloweringRatio, split = ~Li,
            line = list(color='blue',width = 2)) %>% 
  layout(font=list(size=15), title=list(text='(c)',y=0.9, x=0.2),
         scene=list(aspectratio=list(x=1.2,y=1,z=1),
                    xaxis=list(title='<i>L<sub>T</sub></i> (%)', autotick=F,
                               dtick=2.5, range=c(2, 23)),
                    yaxis=list(title='Year'),
                    zaxis=list(title='Flowering synchrony', autotick=F,
                               dtick=0.1, range=c(0.01,0.58)), 
                    camera = list(eye = list(x = 1.3, y = 1.3, z = 1.3),
                                  center = list(x=0, y=0, z=-0.25)))) %>% hide_legend()

Figure 7

Polim_f %>% filter(Year>1983) %>%
  group_by(Rc, Li, Bu, run) %>% summarise(CV=sd(AcornW)/mean(AcornW)) %>%
  group_by(Rc, Li, Bu) %>% summarise(meanCV=mean(CV), seCV=sd(CV)/sqrt(5)) %>%
  arrange(desc(meanCV)) %>%  head(5) %>% left_join(Polim_f) %>%
  group_by(Year, parameter, meanCV) %>%
  summarise(Acorn=mean(AcornW), SD_Acorn=sd(AcornW)) %>%
  convert(fct(parameter)) %>%
  mutate(parameter=fct_reorder(parameter, desc(meanCV))) %>%
  filter(Year>1983) %>%
  group_by(parameter) %>%
  summarize(ratio=SD_Acorn/Acorn, Acorn=(Acorn-mean(Acorn))/sd(Acorn), #ratio to scale later SD
            SD_Acorn=(SD_Acorn-mean(SD_Acorn))/sd(SD_Acorn), #standardized to SD
            SD_Acorn=SD_Acorn*ratio, Year=Year) %>%
  ggplot(aes(x=Year))+
  geom_ribbon(aes(fill=parameter,
                  ymin=Acorn-SD_Acorn,
                  ymax=Acorn+SD_Acorn), alpha=0.1)+
  geom_line(aes(y=Acorn, color=parameter))+
  labs(title='(c)', y='Acorn volume (SD)')+
  theme_bw()+
  theme(legend.title = element_blank(),
        legend.position = c(0.5, 0.88),
        legend.background=element_rect(fill='transparent', size=NULL),
        legend.text.align = 0,
        axis.title.y = element_blank())+
  scale_fill_discrete(labels=c(
    '3_0.175_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=17.5%,',italic('B'[u]),'=4% CV=0.65')),
    '3_0.05_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=5%,',italic('B'[u]),'=4% CV=0.64')),
    '3_0.025_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=2.5%,',italic('B'[u]),'=4% CV=0.63')),
    '3_0.2_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=20%,',italic('B'[u]),'=4% CV=0.62')),
    '3_0.15_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=15%,',italic('B'[u]),'=4% CV=0.61'))))+
  scale_color_discrete(labels=c(
    '3_0.175_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=17.5%,',italic('B'[u]),'=4% CV=0.65')),
    '3_0.05_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=5%,',italic('B'[u]),'=4% CV=0.64')),
    '3_0.025_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=2.5%,',italic('B'[u]),'=4% CV=0.63')),
    '3_0.2_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=20%,',italic('B'[u]),'=4% CV=0.62')),
    '3_0.15_0.04'=
      expression(paste(italic('R'[c]),'=3, ', italic('L'[T]), '=15%,',italic('B'[u]),'=4% CV=0.61'))))+
  guides(fill=guide_legend(ncol=2))+
  geom_line(data=Yacorns, aes(y=(weight-mean(weight, na.rm = T))/sd(weight, na.rm = T)),
            group=1)+
  geom_point(data=Yacorns, aes(y=(weight-mean(weight, na.rm = T))/sd(weight, na.rm = T)))+
  coord_cartesian(ylim=(c(-2.7, 5)))

## Warning: Removed 6 rows containing missing values (geom_point).

Figure 8

#RB+Pollen
Polim_kappa <- 
  Polim_f %>% filter(Year>1983) %>%  select(parameter,Rc, Bu, Li, Year, run, masting) %>% 
  left_join(Yacorns[,c(1,8,9)], by='Year', suffix=c('_sim','_obs')) %>% 
  select(parameter, Rc, Bu, Li,run, Year, masting_sim, mastingW) %>% 
  group_by(parameter,Rc, Bu, Li, run) %>% nest() %>% #by run and parameter
  mutate(data2=map(data, drop_y)) %>% #dropping year from nested frames
  mutate(tables=map(data2,table)) %>% #getting into good table format
  mutate(confKap2=map(tables, fun_Kap2), estimate=map(confKap2, est_Kap),
         p=map(confKap2, p_Kap)) %>% convert(num(estimate,p)) %>% 
  group_by(parameter,Rc, Bu, Li) %>% summarise(est=mean(estimate)) %>% 
  mutate(Li=Li*100, Bu=Bu*100) %>% 
  convert(fct(Bu))

Polim_kappa %>% filter(Bu %in% c('3','4','5')) %>% 
  mutate(Bu=fct_relevel(Bu, "3",'4','5')) %>%
  ggplot(aes(x=jitter(Li, factor = 0.4), y=est, color=Bu)) +
  geom_point()+
  geom_smooth(se=F, size=0.5)+
  labs(title='(c)',x=expression(paste(italic('L'[T]),' (%)')),y='Kappa')+
  scale_color_manual(values=c('black','blue','red'))+
  theme_bw()+
  theme(legend.position = 'none')+
  coord_cartesian(ylim=c(-0.3,0.35))

Figure 9

#RB+Pollen
Polim_group <- Polim_f %>% filter(Year>2004 &Year <2018 &Year!=2010) %>%
  group_by(parameter, Rc, Bu, Li, Year) %>% summarise(AcornW=mean(AcornW)) %>% #mean from 5 runs
  left_join(Yseed, by=c('Year'='year')) %>% 
  mutate(error=(AcornW*1000*2.19)-(weight*0.8)) %>% 
  group_by(parameter, Rc, Bu, Li) %>% 
  summarise(RMSE=sqrt(mean(error^2, na.rm = T)), 
            MAE=mean(abs(error), na.rm=T)) %>% ungroup() %>% 
  mutate(Li=Li*100, Bu=Bu*100) %>% 
  select(Rc:RMSE) %>% convert(fct(Bu))

Polim_group %>% filter(Bu %in% c('3','4','5')) %>% 
  ggplot(aes(x=jitter(Li, factor = 0.4), y=RMSE, color=Bu)) +
  geom_point(size=1)+
  geom_smooth(se=F, size=0.5)+
  labs(title='(c)',x=expression(paste(italic('L'[T]),' (%)')), y='RMSE (kg)')+
  scale_color_manual(values=c('black','blue','red'))+
  theme_bw()+
  theme(legend.position = 'none')+
  coord_cartesian(ylim=c(40,140))

Appendix S4

types <- c("Number" = 'dashed', "Mass" = 'solid')  
Yacorns %>% mutate(CV=AcornSD/AcornMean, CV_W=AcornSDW/AcornMeanW) %>%
  tail(16) %>%
  ggplot(aes(x=Year))+
  geom_line(aes(y=CV, linetype ='Number'))+
  geom_line(aes(y=CV_W, linetype='Mass'))+
  theme_classic()+
  labs(x='End of 20-year window')+
  theme(panel.grid.major.y = element_line(),
        legend.position = 'none',
        axis.text = element_text(size=13),
        axis.title = element_text(size=15))+
  scale_linetype_manual(values = types)+
  scale_x_continuous(breaks=c(2003,2006, 2009, 2012, 2015, 2018))

Appendix S5

## Variability
#RB+Pollen
Polim_f %>% filter(Year>1983) %>% 
  group_by(Rc, Li, Bu, run) %>% summarise(CV=sd(AcornW)/mean(AcornW)) %>% 
  group_by(Rc, Li, Bu) %>% summarise(meanCV=mean(CV), seCV=sd(CV)/sqrt(5)) %>% 
  arrange(desc(meanCV)) %>%  head(5) 

## `summarise()` has grouped output by 'Rc', 'Li', 'Bu'. You can override using the `.groups` argument.

## `summarise()` has grouped output by 'Rc', 'Li'. You can override using the `.groups` argument.

## Simultaneity
#RB+Pollen
Polim_f %>% filter(Year>1983) %>%  select(parameter, Year, run, masting) %>% 
  left_join(Yacorns[,c(1,8,9)], by='Year', suffix=c('_sim','_obs')) %>% 
  select(parameter, run, Year, masting_sim, mastingW) %>% 
  group_by(parameter, run) %>% nest() %>% #by run and parameter
  mutate(data2=map(data, drop_y)) %>% #dropping year from nested frames
  mutate(tables=map(data2,table)) %>% #getting into good table format
  mutate(confKap2=map(tables, fun_Kap2), estimate=map(confKap2, est_Kap),
         p=map(confKap2, p_Kap)) %>% convert(num(estimate,p)) %>% 
  group_by(parameter) %>% summarise(est=mean(estimate)) %>%  
  arrange(desc(est)) %>% head(5)

## Effect size
Yacorn_freq <- Yacorns %>% filter(!is.na(freq)) %>% 
  select(Year, weight, freqW)

#RB+Pollen
Polim_f %>% filter(Year>2002) %>% #end of 20 year running window
  left_join(Yacorn_freq) %>% 
  group_by(parameter, run) %>% nest() %>%
  mutate(wilTest=map(data, wil_test_eff)) %>% 
  unnest(wilTest) %>% group_by(parameter) %>% summarise(eff=mean(effsize)) %>% 
  arrange(eff) %>% head(5)

## Joining, by = "Year"

## RMSE
#RB+Pollen
Polim_f %>% filter(Year>2004 &Year <2018 &Year!=2010) %>%
  group_by(parameter, Year) %>% summarise(AcornW=mean(AcornW)) %>% #mean from 5 runs
  left_join(Yseed, by=c('Year'='year')) %>% 
  mutate(error=(AcornW*1000*2.19)-(weight*0.8)) %>% 
  group_by(parameter) %>% 
  summarise(RMSE=sqrt(mean(error^2, na.rm = T))) %>% 
  arrange(RMSE) %>% head(5)

## `summarise()` has grouped output by 'parameter'. You can override using the `.groups` argument.